#unweighted rep on data set three

library(Amelia)
library(MatchIt)
library(cem)
library(foreign)
library(survey)
library("Zelig")

#import

data <- read.csv("three_final3.csv", na.string=c(" "))

#original replication
#define survey weights, strata and cluster

weights <- svydesign(id=~CPSUM, strata=~CSTRATM, weights=~PATWT, data=data)

#reproduce table 1 - frequencies
output.t.f <- svyby(~any_image+adv_image+MRI+CATSCAN, ~female, weights, svymean)
rownames(output.t.f) <- c("Male","Female")
output.t.a <-svyby(~any_image+adv_image+MRI+CATSCAN, ~age1, weights, svymean)
rownames(output.t.a) <- c("< 18","18-45", "46-64", ">64")
output.t.bl <-svyby(~any_image+adv_image+MRI+CATSCAN, ~black, weights, svymean)
rownames(output.t.bl) <- c("Nonblack", "Black")
output.t.ins <-svyby(~any_image+adv_image+MRI+CATSCAN, ~insurance_cat, weights, svymean, na.rm=TRUE)
rownames(output.t.ins) <- c("Private", "Medicare","Medicaid", "Other", "None")
output.t.md_cat <-svyby(~any_image+adv_image+MRI+CATSCAN, ~mdtype_cat, weights, svymean)
rownames(output.t.md_cat) <- c("Surgical","Primary","Medical")
output.t.pract_set <-svyby(~any_image+adv_image+MRI+CATSCAN, ~prac_setting, weights, svymean)
rownames(output.t.pract_set) <- c("Private Office","Community Health Center","HMO", "Free-standing Clinic","Other")
output.t.md_own <-svyby(~any_image+adv_image+MRI+CATSCAN, ~md_owns, weights, svymean, na.rm=TRUE)
rownames(output.t.md_own) <- c("No","Yes") 
output.t.imgr <-svyby(~any_image+adv_image+MRI+CATSCAN, ~image_result, weights, svymean, na.rm=TRUE)
rownames(output.t.imgr) <- c("No","Yes")
output.t.imgv <-svyby(~any_image+adv_image+MRI+CATSCAN, ~image_view, weights, svymean, na.rm=TRUE)
rownames(output.t.imgv) <- c("No","Yes")
#note this produces warning messages that do not affect the results

#reproduce table 1 - n's
n.output.t.f <- svyby(~any_image+adv_image+MRI+CATSCAN, ~female, weights, unwtd.count)
rownames(n.output.t.f) <- c("Male","Female")
n.output.t.a <-svyby(~any_image+adv_image+MRI+CATSCAN, ~age1, weights, unwtd.count, na.rm=TRUE)
rownames(n.output.t.a) <- c("< 18","18-45", "46-64", ">64")
n.output.t.bl <-svyby(~any_image+adv_image+MRI+CATSCAN, ~black, weights, unwtd.count)
rownames(n.output.t.bl) <- c("Nonblack", "Black")
n.output.t.ins <-svyby(~any_image+adv_image+MRI+CATSCAN, ~insurance_cat, weights, unwtd.count, na.rm=TRUE)
rownames(n.output.t.ins) <- c("Private", "Medicare","Medicaid", "Other", "None")
n.output.t.md_cat <-svyby(~any_image+adv_image+MRI+CATSCAN, ~mdtype_cat, weights, unwtd.count)
rownames(n.output.t.md_cat) <- c("Surgical","Primary","Medical")
n.output.t.pract_set <-svyby(~any_image+adv_image+MRI+CATSCAN, ~prac_setting, weights, unwtd.count)
rownames(n.output.t.pract_set) <- c("Private Office","Community Health Center","HMO", "Free-standing Clinic","Other")
n.output.t.md_own <-svyby(~any_image+adv_image+MRI+CATSCAN, ~md_owns, weights, unwtd.count, na.rm=TRUE)
rownames(n.output.t.md_own) <- c("No","Yes") 
n.output.t.imgr <-svyby(~any_image+adv_image+MRI+CATSCAN, ~image_result, weights, unwtd.count, na.rm=TRUE)
rownames(n.output.t.imgr) <- c("No","Yes")
n.output.t.imgv <-svyby(~any_image+adv_image+MRI+CATSCAN, ~image_view, weights, unwtd.count, na.rm=TRUE)
rownames(n.output.t.imgv) <- c("No","Yes")

#list both frequencies and n's
list(round(output.t.f[,2:5], digits=3)*100, round(output.t.a[,2:5], digits=3)*100,round(output.t.bl[,2:5], digits=3)*100, round(output.t.ins[1:5,2:5], digits=3)*100,round(output.t.md_cat[,2:5], digits=3)*100,round(output.t.pract_set[,2:5], digits=3)*100,round(output.t.md_own[1:2,2:5], digits=3)*100,round(output.t.imgr[1:2,2:5], digits=3)*100,round(output.t.imgv[1:2,2:5],digits=3)*100)
list(n.output.t.f[,1:2],n.output.t.a[,1:2],n.output.t.bl[,1:2], n.output.t.ins[1:2], n.output.t.md_cat[,1:2][,1:2], n.output.t.pract_set[,1:2], n.output.t.md_own[,1:2], n.output.t.imgr[,1:2], n.output.t.imgv[,1:2])

##tables 2 and 3
# replicate odds ratios

model.1.r <- zelig(any_image~female+age_lt45+age_lt64+age_gt64+black+hispanic+poverty_cat+urban_cat+seen+ins_medicare+ins_medicaid+ins_other+ins_none+md_primary+md_medical+setting_ctty+setting_hmo+setting_free+setting_other+md_owns+solo_cat+prepaid+hosp_owns+cost_profiling+image_result, model="logit.survey", data=data, ids=data$CPSUM, strata=data$CSTRATM)
odds.1.r <- exp(model.1.r$coefficients)
odds.1.r <- round(odds.1.r, digits=2)
odds.1.r <- as.matrix(odds.1.r)
no.image.1.r <- setx(model.1.r, fn=list(numeric=median), image_result=0)
image.1.r <- setx(model.1.r, fn=list(numeric=median), age_lt64=1,image_result=1)
image.1.r.sim <- sim(model.1.r, x=no.image.1.r, x1=image.1.r)
summary(image.1.r.sim)
plot(image.1.r.sim)

model.1.v <- zelig(any_image~female+age_lt45+age_lt64+age_gt64+black+hispanic+poverty_cat+urban_cat+seen+ins_medicare+ins_medicaid+ins_other+ins_none+md_primary+md_medical+setting_ctty+setting_hmo+setting_free+setting_other+md_owns+solo_cat+prepaid+hosp_owns+cost_profiling+image_view, model="logit.survey", data=data, weights=data$PATWT, ids=data$CPSUM, strata=data$CSTRATM)
odds.1.v <- exp(model.1.v$coefficients)
odds.1.v <-round(odds.1.v, digits=2)
odds.1.v <- as.matrix(odds.1.v)

no.image.1.v <- setx(model.1.v, fn=list(numeric=median), age_lt64=1, image_view=0)
image.1.v <- setx(model.1.v, fn=list(numeric=median), age_lt64=1, image_view=1)
image.1.v.sim <- sim(model.1.v, x=no.image.1.v, x1=image.1.v)
summary(image.1.v.sim)
plot(image.1.v.sim)

model.2.r <- zelig(adv_image~female+age_lt45+age_lt64+age_gt64+black+hispanic+poverty_cat+urban_cat+seen+ins_medicare+ins_medicaid+ins_other+ins_none+md_primary+md_medical+setting_ctty+setting_hmo+setting_free+setting_other+md_owns+solo_cat+prepaid+hosp_owns+cost_profiling+image_result, model="logit.survey", data=data, weights=data$PATWT, ids=data$CPSUM, strata=data$CPSTRATM)
odds.2.r <- exp(model.2.r$coefficients)
odds.2.r <-round(odds.2.r, digits=2)
odds.2.r <- as.matrix(odds.2.r)

no.image.2.r <- setx(model.2.r, fn=list(numeric=median),age_lt64=1, image_result=0)
image.2.r <- setx(model.2.r, fn=list(numeric=median), age_lt64=1,image_result=1)
image.2.r.sim <- sim(model.2.r, x=no.image.2.r, x1=image.2.r)
summary(image.2.r.sim)
plot(image.2.r.sim)

model.2.v <- zelig(adv_image~female+age_lt45+age_lt64+age_gt64+black+hispanic+poverty_cat+urban_cat+seen+ins_medicare+ins_medicaid+ins_other+ins_none+md_primary+md_medical+setting_ctty+setting_hmo+setting_free+setting_other+md_owns+solo_cat+prepaid+hosp_owns+cost_profiling+image_view, model="logit.survey", data=data, weights=data$PATWT, ids=data$CPSUM, strata=data$CPSTRATM)
odds.2.v <- exp(model.2.v$coefficients)
odds.2.v <- round(odds.2.v, digit=2)
odds.2.v <- as.matrix(odds.2.v)

no.image.2.v <- setx(model.2.v, fn=list(numeric=median), age_lt64=1,image_view=0)
image.2.v <- setx(model.2.v, fn=list(numeric=median), age_lt64=1,image_view=1)
image.2.v.sim <- sim(model.2.v, x=no.image.2.v, x1=image.2.v)
summary(image.2.v.sim)
plot(image.2.v.sim)

odds.all <- cbind(odds.1.r, odds.1.v, odds.2.r, odds.2.v)
colnames(odds.all) <- c("Any image-Results", "Any image-Views", "Adv. image-Results", "Adv. image-Views")
odds.all

#list model summary to show significance based on t-tests
summary(model.1.r)
summary(model.1.v)
summary(model.2.r)
summary(model.2.v)

#check original imbalance

data.vars <- subset(data, select=c("any_image","adv_image","female","age_lt45","age_lt64","age_gt64","black","hispanic","poverty_cat","urban_cat","seen","ins_medicare","ins_medicaid","ins_other","ins_none","md_primary","md_medical","setting_ctty","setting_hmo","setting_free","setting_other","md_owns","solo_cat","prepaid","hosp_owns","cost_profiling","image_result", "image_view", "ASTHMA", "CANCER", "CEBVD", "CHF", "COPD", "DEPRN", "DIABETES", "HYPLIPID", "HTN", "IHD", "OSTPRSIS", "OBESITY", "PT", "ORTHO", "PCCQOC","dmed_1","dmed_2","dmed_3","dmed_4","ccchron_con_1","ccchron_con_2","ccchron_con_3","ccchron_con_4","rdregion_2","rdregion_3","rdregion_4","dmmajor_2","dmmajor_3","dmmajor_4","dmmajor_5"))
imbalance(group=data.vars$image_result, data=data.vars, drop=c("any_image","adv_image","image_result", "image_view"))

#checking a few frequencies
attach(data)
female.table <- table(female, image_result)
black.talbe <- table(black, image_result)
insnone.table <-table(ins_none, image_result)
solo.table <-table(solo_cat, image_result)
prop.table(solo.table,2)

#MI
set.seed(123)
data.imp <- amelia(data.vars, m=5, frontend=FALSE, ords=c("any_image","adv_image","female","age_lt45","age_lt64","age_gt64","black","hispanic","poverty_cat","urban_cat","seen","ins_medicare","ins_medicaid","ins_other","ins_none","md_primary","md_medical","setting_ctty","setting_hmo","setting_free","setting_other","md_owns","solo_cat","prepaid","hosp_owns","cost_profiling","image_result", "image_view", "ASTHMA", "CANCER", "CEBVD", "CHF", "COPD", "DEPRN", "DIABETES", "HYPLIPID", "HTN", "IHD", "OSTPRSIS", "OBESITY", "PT", "ORTHO", "PCCQOC","dmed_1","dmed_2","dmed_3","dmed_4","ccchron_con_1","ccchron_con_2","ccchron_con_3","ccchron_con_4","rdregion_2","rdregion_3","rdregion_4","dmmajor_2","dmmajor_3","dmmajor_4","dmmajor_5"))
compare.density(data.imp, var="image_result")
summary(data.imp)
imbalance(group=data.imp$imputations[1]$imp1$image_view, data=data.imp$imputations[1]$imp1,drop=c("any_image","adv_image","image_result", "image_view") )
imbalance(group=data.imp$imputations[2]$imp2$image_view, data=data.imp$imputations[1]$imp1,drop=c("any_image","adv_image","image_result", "image_view") )
imbalance(group=data.imp$imputations[3]$imp3$image_view, data=data.imp$imputations[1]$imp1,drop=c("any_image","adv_image","image_result", "image_view") )
imbalance(group=data.imp$imputations[4]$imp4$image_view, data=data.imp$imputations[1]$imp1,drop=c("any_image","adv_image","image_result", "image_view") )
imbalance(group=data.imp$imputations[5]$imp5$image_view, data=data.imp$imputations[1]$imp1,drop=c("any_image","adv_image","image_result", "image_view") )

#matching with image result as treatment
dataframe.imp1.1 <- as.data.frame(data.imp$imputations[1])
dataframe.imp1.2 <- as.data.frame(data.imp$imputations[2])
dataframe.imp1.3 <- as.data.frame(data.imp$imputations[3])
dataframe.imp1.4 <- as.data.frame(data.imp$imputations[4])
dataframe.imp1.5 <- as.data.frame(data.imp$imputations[5])

data.cem.match.imp1.1 <- matchit(imp1.image_result~imp1.md_primary+imp1.md_medical+imp1.setting_ctty+imp1.setting_hmo+imp1.setting_free+imp1.setting_other+imp1.md_owns+imp1.solo_cat+imp1.prepaid+imp1.hosp_owns+imp1.cost_profiling+imp1.PCCQOC, data=dataframe.imp1.1, method="cem", drop =c("any_image","adv_image","female","age_lt45","age_lt64","age_gt64","black","hispanic","poverty_cat","urban_cat","ins_medicare","ins_medicaid","ins_other","ins_none","image_result", "image_view", "ASTHMA", "CANCER", "CEBVD", "CHF", "COPD", "DEPRN", "DIABETES", "HYPLIPID", "HTN", "IHD", "OSTPRSIS", "OBESITY", "PT", "ORTHO","dmed_1","dmed_2","dmed_3","dmed_4","ccchron_con_1","ccchron_con_2","ccchron_con_3","ccchron_con_4","rdregion_2","rdregion_3","rdregion_4","dmmajor_2","dmmajor_3","dmmajor_4","dmmajor_5"))
dataframe.cem.match.imp1.1 <- match.data(data.cem.match.imp1.1)
data.cem.match.imp1.2 <- matchit(imp2.image_result~imp2.md_primary+imp2.md_medical+imp2.setting_ctty+imp2.setting_hmo+imp2.setting_free+imp2.setting_other+imp2.md_owns+imp2.solo_cat+imp2.prepaid+imp2.hosp_owns+imp2.cost_profiling+imp2.PCCQOC, data=dataframe.imp1.2, method="cem", drop =c("any_image","adv_image","female","age_lt45","age_lt64","age_gt64","black","hispanic","poverty_cat","urban_cat","ins_medicare","ins_medicaid","ins_other","ins_none","image_result", "image_view", "ASTHMA", "CANCER", "CEBVD", "CHF", "COPD", "DEPRN", "DIABETES", "HYPLIPID", "HTN", "IHD", "OSTPRSIS", "OBESITY", "chron_conditions", "MED", "REGION", "PT", "ORTHO", "MAJOR"))
dataframe.cem.match.imp1.2 <- match.data(data.cem.match.imp1.2)
data.cem.match.imp1.3 <- matchit(imp3.image_result~imp3.md_primary+imp3.md_medical+imp3.setting_ctty+imp3.setting_hmo+imp3.setting_free+imp3.setting_other+imp3.md_owns+imp3.solo_cat+imp3.prepaid+imp3.hosp_owns+imp3.cost_profiling+imp3.PCCQOC, data=dataframe.imp1.3, method="cem", drop =c("any_image","adv_image","female","age_lt45","age_lt64","age_gt64","black","hispanic","poverty_cat","urban_cat","ins_medicare","ins_medicaid","ins_other","ins_none","image_result", "image_view", "ASTHMA", "CANCER", "CEBVD", "CHF", "COPD", "DEPRN", "DIABETES", "HYPLIPID", "HTN", "IHD", "OSTPRSIS", "OBESITY", "chron_conditions", "MED", "REGION", "PT", "ORTHO", "MAJOR"))
dataframe.cem.match.imp1.3 <- match.data(data.cem.match.imp1.3)
data.cem.match.imp1.4 <- matchit(imp4.image_result~imp4.md_primary+imp4.md_medical+imp4.setting_ctty+imp4.setting_hmo+imp4.setting_free+imp4.setting_other+imp4.md_owns+imp4.solo_cat+imp4.prepaid+imp4.hosp_owns+imp4.cost_profiling+imp4.PCCQOC, data=dataframe.imp1.4, method="cem", drop =c("any_image","adv_image","female","age_lt45","age_lt64","age_gt64","black","hispanic","poverty_cat","urban_cat","ins_medicare","ins_medicaid","ins_other","ins_none","image_result", "image_view", "ASTHMA", "CANCER", "CEBVD", "CHF", "COPD", "DEPRN", "DIABETES", "HYPLIPID", "HTN", "IHD", "OSTPRSIS", "OBESITY", "chron_conditions", "MED", "REGION", "PT", "ORTHO", "MAJOR"))
dataframe.cem.match.imp1.4 <- match.data(data.cem.match.imp1.4)
data.cem.match.imp1.5 <- matchit(imp5.image_result~imp5.md_primary+imp5.md_medical+imp5.setting_ctty+imp5.setting_hmo+imp5.setting_free+imp5.setting_other+imp5.md_owns+imp5.solo_cat+imp5.prepaid+imp5.hosp_owns+imp5.cost_profiling+imp5.PCCQOC, data=dataframe.imp1.5, method="cem", drop =c("any_image","adv_image","female","age_lt45","age_lt64","age_gt64","black","hispanic","poverty_cat","urban_cat","ins_medicare","ins_medicaid","ins_other","ins_none","image_result", "image_view", "ASTHMA", "CANCER", "CEBVD", "CHF", "COPD", "DEPRN", "DIABETES", "HYPLIPID", "HTN", "IHD", "OSTPRSIS", "OBESITY", "chron_conditions", "MED", "REGION", "PT", "ORTHO", "MAJOR"))
dataframe.cem.match.imp1.5 <- match.data(data.cem.match.imp1.5)

summary(data.cem.match.imp1.1)
summary(data.cem.match.imp1.2)
summary(data.cem.match.imp1.3)
summary(data.cem.match.imp1.4)
summary(data.cem.match.imp1.5)

#create frequency tables for matched/imputed data

stderr <- function(x) sqrt(var(x)/length(x))

mimean1 <-aggregate(dataframe.cem.match.imp1.1, list(Image_result=dataframe.cem.match.imp1.1$imp1.image_result), mean)
mise1 <- aggregate(dataframe.cem.match.imp1.1, list(Image_result=dataframe.cem.match.imp1.1$imp1.image_result), stderr)
mimean2 <- aggregate(dataframe.cem.match.imp1.2, list(Image_result=dataframe.cem.match.imp1.2$imp2.image_result), mean)
mise2 <-aggregate(dataframe.cem.match.imp1.2, list(Image_result=dataframe.cem.match.imp1.2$imp2.image_result), stderr)
mimean3 <-aggregate(dataframe.cem.match.imp1.3, list(Image_result=dataframe.cem.match.imp1.3$imp3.image_result), mean)
mise3 <-aggregate(dataframe.cem.match.imp1.3, list(Image_result=dataframe.cem.match.imp1.3$imp3.image_result), stderr)
mimean4 <-aggregate(dataframe.cem.match.imp1.4, list(Image_result=dataframe.cem.match.imp1.4$imp4.image_result), mean)
mise4 <-aggregate(dataframe.cem.match.imp1.4, list(Image_result=dataframe.cem.match.imp1.4$imp4.image_result), stderr)
mimean5 <-aggregate(dataframe.cem.match.imp1.5, list(Image_result=dataframe.cem.match.imp1.5$imp5.image_result), mean)
mise5 <-aggregate(dataframe.cem.match.imp1.5, list(Image_result=dataframe.cem.match.imp1.5$imp5.image_result), stderr)

colnames(mimean2) <- colnames(mimean1)
colnames(mimean3) <- colnames(mimean1)
colnames(mimean4) <- colnames(mimean1)
colnames(mimean5) <- colnames(mimean1)

colnames(mise2) <- colnames(mimean1)
colnames(mise3) <- colnames(mimean1)
colnames(mise4) <- colnames(mimean1)
colnames(mise5) <- colnames(mimean1)

mimeans0 <- rbind(mimean1[1,],mimean2[1,],mimean3[1,],mimean4[1,],mimean5[1,])
mimeans1 <- rbind(mimean1[2,],mimean2[2,],mimean3[2,],mimean4[2,],mimean5[2,])
mises0 <- rbind(mise1[1,],mise2[1,],mise3[1,],mise4[1,],mise5[1,])
mises1 <- rbind(mise1[2,],mise2[2,],mise3[2,],mise4[2,],mise5[2,])

no.img.frq <- mi.meld(mimeans0, mises0)
yes.img.frq <- mi.meld(mimeans1, mises1)
cbind(round(t(no.img.frq$q.mi), digits=4),round(t(yes.img.frq$q.mi), digits=4))

#check balance of matched and imputed data

data.cem.imp.1$match1$imbalance
data.cem.imp.1$match2$imbalance
data.cem.imp.1$match3$imbalance
data.cem.imp.1$match4$imbalance
data.cem.imp.1$match5$imbalance

#####comparing restricted and unrestricted models based on matched data 

imputed.matched.reg1.adj <- zelig(imp1.any_image ~ imp1.image_result
, model="logit", data=dataframe.cem.match.imp1.1)

print(summary(imputed.matched.reg1.adj), digits = 2)

loglik <- -1/2*deviance(imputed.matched.reg1.adj)

loglik

imputed.matched.reg1.adj.pretxt <- zelig(imp1.any_image ~ imp1.image_result + imp1.seen + imp1.md_primary + imp1.md_medical + imp1.setting_ctty + 
imp1.setting_hmo + imp1.setting_free + imp1.setting_other + imp1.md_owns + imp1.solo_cat + imp1.prepaid + imp1.hosp_owns + imp1.cost_profiling + imp1.PCCQOC
, data=dataframe.cem.match.imp1.1, model="logit")
summary(imputed.matched.reg1.adj.pretxt)

print(summary(imputed.matched.reg1.adj.pretxt), digits = 2)

loglik.pretxt <- -1/2*deviance(imputed.matched.reg1.adj.pretxt)

loglik.pretxt

#likelihood ratio test
lrt <- 2*(loglik.pretxt - loglik)
1 - pchisq(lrt, df = 13)

imputed.matched.reg1.adj.all <- zelig(imp1.any_image ~ imp1.image_result + imp1.seen + imp1.md_primary + imp1.md_medical + imp1.setting_ctty + 
imp1.setting_hmo + imp1.setting_free + imp1.setting_other + imp1.md_owns + imp1.solo_cat + imp1.prepaid + imp1.hosp_owns + imp1.cost_profiling + imp1.PCCQOC 
+ imp1.REGION + imp1.female + imp1.age_lt45 + imp1.age_lt64 + imp1.age_gt64 + imp1.black + imp1.hispanic + imp1.poverty_cat + imp1.urban_cat + imp1.ins_medicare 
+ imp1.ins_medicaid + imp1.ins_other + imp1.ins_none + imp1.ASTHMA + imp1.CANCER + imp1.CEBVD + imp1.CHF + imp1.COPD + imp1.DEPRN + imp1.DIABETES + 
imp1.HYPLIPID + imp1.HTN + imp1.IHD + imp1.OSTPRSIS + imp1.OBESITY + imp1.MED + imp1.REGION + imp1.PT + imp1.ORTHO + imp1.MAJOR
, data=dataframe.cem.match.imp1.1, model="logit")

summary(imputed.matched.reg1.adj.all)

print(summary(imputed.matched.reg1.adj.all), digits = 2)

loglik.all <- -1/2*deviance(imputed.matched.reg1.adj.all)

loglik.all

#likelihood ratio test

lrt.all.pretxt <- 2*(loglik.all - loglik.pretxt)
1 - pchisq(lrt.all.pretxt, df = 30)

#evaluating model fit using ROC curve

rocplot(imputed.matched.reg1.adj$y, imputed.matched.reg1.adj.pretxt$y, fitted(imputed.matched.reg1.adj), fitted(imputed.matched.reg1.adj.pretxt))

rocplot(imputed.matched.reg1.adj.all$y, imputed.matched.reg1.adj.pretxt$y, fitted(imputed.matched.reg1.adj.all), fitted(imputed.matched.reg1.adj.pretxt))


#QOI - evaluating first differences for above models

media <- dataframe.cem.match.imp1.1[dataframe.cem.match.imp1.1$imp1.image_result == 1,]

median.cov <- media[,c("imp1.image_result", "imp1.seen", "imp1.md_primary", "imp1.md_medical", "imp1.setting_ctty", 
"imp1.setting_hmo", "imp1.setting_free", "imp1.setting_other", "imp1.md_owns", "imp1.solo_cat", "imp1.prepaid", "imp1.hosp_owns", 
"imp1.cost_profiling", "imp1.PCCQOC", "imp1.REGION", "imp1.female", "imp1.age_lt45", "imp1.age_lt64", "imp1.age_gt64", "imp1.black", 
"imp1.hispanic", "imp1.poverty_cat", "imp1.urban_cat", "imp1.ins_medicare", "imp1.ins_medicaid", "imp1.ins_other", "imp1.ins_none", 
"imp1.ASTHMA", "imp1.CANCER", "imp1.CEBVD", "imp1.CHF", "imp1.COPD", "imp1.DEPRN",
"imp1.DIABETES", "imp1.HYPLIPID", "imp1.HTN", "imp1.IHD", "imp1.OSTPRSIS", "imp1.OBESITY", "imp1.MED", "imp1.REGION", "imp1.PT", 
"imp1.ORTHO", "imp1.MAJOR")]

med.ian <- apply(median.cov, 2, median)
med.ian

#1st difference - unadjusted model

image_result.yes <- setx(imputed.matched.reg1.adj, imp1.image_result = 1)


image_result.no <- setx(imputed.matched.reg1.adj, imp1.image_result = 0)

s.out <- sim(imputed.matched.reg1.adj, x = image_result.no, x1 = image_result.yes)

summary(s.out)

plot(s.out)

#1st difference - pretxt var adjusted model

image_result.pre.yes <- setx(imputed.matched.reg1.adj.pretxt, imp1.image_result = 1, imp1.seen=1, imp1.md_primary = 0,
imp1.md_medical = 0, imp1.setting_ctty = 0, imp1.setting_hmo = 0,  imp1.setting_free = 0, 
imp1.setting_other = 0, imp1.md_owns = 1, imp1.solo_cat = 0, imp1.prepaid = 0, imp1.hosp_owns = 0,
imp1.cost_profiling = 0, imp1.PCCQOC = 2)


image_result.pre.no <- setx(imputed.matched.reg1.adj.pretxt, imp1.image_result = 0, imp1.seen=1, imp1.md_primary = 0,
imp1.md_medical = 0, imp1.setting_ctty = 0, imp1.setting_hmo = 0,  imp1.setting_free = 0, 
imp1.setting_other = 0, imp1.md_owns = 1, imp1.solo_cat = 0, imp1.prepaid = 0, imp1.hosp_owns = 0,
imp1.cost_profiling = 0, imp1.PCCQOC = 2)

s.out.pre <- sim(imputed.matched.reg1.adj.pretxt, x = image_result.pre.no, x1 = image_result.pre.yes)

summary(s.out.pre)

plot(s.out.pre)

#1st difference - all var adjusted model

image_result.all.yes <- setx(imputed.matched.reg1.adj.all, imp1.image_result = 1, imp1.seen=1, imp1.md_primary = 0,
imp1.md_medical = 0, imp1.setting_ctty = 0, imp1.setting_hmo = 0,  imp1.setting_free = 0, 
imp1.setting_other = 0, imp1.md_owns = 1, imp1.solo_cat = 0, imp1.prepaid = 0, imp1.hosp_owns = 0,
imp1.cost_profiling = 0, imp1.PCCQOC = 2,
imp1.REGION = 3, imp1.female =1, imp1.age_lt45 = 0, imp1.age_lt64 = 0, imp1.age_gt64 = 0, 
imp1.black = 0, imp1.hispanic = 0, imp1.poverty_cat = 0, imp1.urban_cat = 1, imp1.ins_medicare = 0,
imp1.ins_medicaid = 0, imp1.ins_other = 0, imp1.ins_none = 0, imp1.ASTHMA = 0, imp1.CANCER = 0, 
imp1.CEBVD = 0, imp1.CHF = 0, imp1.COPD = 0, imp1.DEPRN = 0, imp1.DIABETES = 0,  
imp1.HYPLIPID = 0, imp1.HTN = 0, imp1.IHD = 0, imp1.OSTPRSIS = 0, imp1.OBESITY = 0,
imp1.MED = 2, imp1.PT = 0, imp1.ORTHO = 1, imp1.MAJOR = 2)


image_result.all.no <- setx(imputed.matched.reg1.adj.all, imp1.image_result = 0, imp1.seen=1, imp1.md_primary = 0,
imp1.md_medical = 0, imp1.setting_ctty = 0, imp1.setting_hmo = 0,  imp1.setting_free = 0, 
imp1.setting_other = 0, imp1.md_owns = 1, imp1.solo_cat = 0, imp1.prepaid = 0, imp1.hosp_owns = 0,
imp1.cost_profiling = 0, imp1.PCCQOC = 2,
imp1.REGION = 3, imp1.female =1, imp1.age_lt45 = 0, imp1.age_lt64 = 0, imp1.age_gt64 = 0, 
imp1.black = 0, imp1.hispanic = 0, imp1.poverty_cat = 0, imp1.urban_cat = 1, imp1.ins_medicare = 0,
imp1.ins_medicaid = 0, imp1.ins_other = 0, imp1.ins_none = 0, imp1.ASTHMA = 0, imp1.CANCER = 0, 
imp1.CEBVD = 0, imp1.CHF = 0, imp1.COPD = 0, imp1.DEPRN = 0, imp1.DIABETES = 0,  
imp1.HYPLIPID = 0, imp1.HTN = 0, imp1.IHD = 0, imp1.OSTPRSIS = 0, imp1.OBESITY = 0,
imp1.MED = 2, imp1.PT = 0, imp1.ORTHO = 1, imp1.MAJOR = 2)

s.out.all <- sim(imputed.matched.reg1.adj.all, x = image_result.all.no, x1 = image_result.all.yes)

summary(s.out.all)

plot(s.out.all)






